# Set of climatic variables# =========================clim_var <-readRDS(here("data/ClimaticData/NamesSelectedVariables.rds"))# Climatic data# =============# we scale the past and future climatic data at the location of the populations# with the parameters (mean and variance) of the past climatic data.source(here("scripts/functions/generate_scaled_clim_datasets.R"))clim_dfs <-generate_scaled_clim_datasets(clim_var,clim_ref_adj =FALSE)
Adaptive landscape
Projecting adaptive gradient(s) across space
Adaptively enriched genetic space
Following Capblancq and Forester (2021) (and the associated Github repository), for the different sets of SNPs (one set of control SNPs and two sets of putative adaptive loci), a RDA was performed with the set of SNPs as multivariate response and the set of climatic variables as explanatory variables. For the sets of candidate SNPs, the resulting RDA spaces can be considered as ‘adaptively enriched genetic spaces’ (Capblancq and Forester 2021; Capblancq et al. 2023).
Code
runRDA <-function(snp_set, clim_var, clim_df, geno){# Keep the genomic data of the selected SNPsgeno <- geno[,snp_set]# RDA formula form_rda <-paste("geno ~ ", paste(clim_var,collapse=" + "))# run RDA with only the selected SNPsrda_outliers <-rda(as.formula(form_rda), clim_df)}
rda_biplots <-lapply(snp_sets, make_RDAbiplot)grid_RDAbiplots <-plot_grid(rda_biplots[[1]],rda_biplots[[2]],rda_biplots[[3]], nrow=1)ggsave(here("figs/RDA/RDAbiplots_SNPsets.pdf"), grid_RDAbiplots, width =14, height=6)# We save a plot for the Supplementary Information# ================================================plot_grid(rda_biplots[[1]],rda_biplots[[3]],nrow=1) %>%# we keep only the control and candidate SNPsggsave(here("figs/RDA/RDAbiplots_SNPsets_SI.pdf"), ., width =9, height=6)grid_RDAbiplots
Adaptive index across the landscape
We used the scores of the climatic variables along the RDA axes to calculate a genetic-based index for each pixel of the landscape. For each RDA axis, the index is calculated as follows:
\[\sum_{i=1}^{n}a_ib_i\]
where \(i\) is one of the \(n\) climatic variables used in the RDA model, \(a\) is the score of the climatic variable \(i\) along the RDA axis and \(b\) is the standardized value for the climatic variable \(i\) at the focal pixel.
When calculating based on the set of candidate SNPs, this index can be considered as an adaptive index that provides an estimate of the adaptive genetic similarity or difference of all pixels on the landscape as a function of the values of the climatic predictors at that location.
We project the index on maps to visualize the different gradients across the maritime pine range (which can be considered as adaptive gradients when based on the sets of putatively adaptive loci).
Code
# countries borders for the mapworld <-ne_countries(scale ="medium", returnclass ="sf")ai_maps <-lapply(snp_sets, function(snp_set){RDA_proj <- ai_proj[[snp_set$set_code]]RDA_proj <-lapply(RDA_proj, function(x) rasterToPoints(x))# The adaptive index is scaled between 0 and 1for(i in1:length(RDA_proj)){ RDA_proj[[i]][,3] <- (RDA_proj[[i]][,3]-min(RDA_proj[[i]][,3]))/(max(RDA_proj[[i]][,3])-min(RDA_proj[[i]][,3]))}# formatting the dataframes for ggplotTAB_RDA <-as.data.frame(do.call(rbind, RDA_proj[1:2]))colnames(TAB_RDA)[3] <-"value"TAB_RDA$variable <-factor(c(rep("RDA1", nrow(RDA_proj[[1]])), rep("RDA2", nrow(RDA_proj[[2]]))), levels =c("RDA1","RDA2"))# map optionspoint_size=2x_limits =c(-10, 15)y_limits =c(31, 52)ggtitle <- snp_set$set_nameggplot(data = TAB_RDA) +geom_sf(data = world, fill="gray98") +scale_x_continuous(limits = x_limits) +scale_y_continuous(limits = y_limits) +geom_raster(aes(x = x, y = y, fill =cut(value, breaks=seq(0, 1, length.out=11), include.lowest = T))) +scale_fill_viridis_d(alpha =0.8, direction =-1, option ="A") +xlab("Longitude") +ylab("Latitude") +guides(fill=guide_legend(title="Genetic index")) +facet_grid(~ variable) +ggtitle(ggtitle) +theme_bw(base_size =11) +theme(panel.grid =element_blank(), plot.background =element_blank(), panel.background =element_blank(), strip.text =element_text(size=11))})# We save the maps for the Github repository# ==========================================pdf(here("figs/RDA/AdaptiveIndex_maps.pdf"), width=10,height=6)lapply(ai_maps, function(x) x)dev.off()# We save the maps for the Supplementary Information# ==================================================# Candidate SNPspdf(here("figs/RDA/AdaptiveIndexMap_CandidateSNPs_SI.pdf"), width=9,height=4.9)ai_maps[[1]] +theme(plot.title =element_blank())dev.off()# Control SNPspdf(here("figs/RDA/AdaptiveIndexMap_ControlSNPs_SI.pdf"), width=9,height=4.9)ai_maps[[3]] +theme(plot.title =element_blank())dev.off()# Show AI maps# ============lapply(ai_maps, function(x) x)
Genomic offset predictions
Using spatial points
Predicting the genomic offset of the studied populations under future climates.
Code
# Function to calculate the RDA genetic offset for spatial points# ===============================================================# clim_ref = climatic data used to fit the RDA model (the climate-of-origin of the populations)# clim_pred = climatic data used for predictions (so either the future climate of the populations, # or climate of the common gardens or the NFI plots)# weights = if NULL, the adaptive index is not weighted by the relative importance of the RDA axes in # explaining the variance.# K = number of RDA axes used to calculate the genomic offset# snp_set = list with at least the RDA models pred_GO_spatial_points <-function(snp_set, K, clim_ref, clim_pred, clim_var,weights =NULL, CG=F){# weights based on the variance explained by the different axesweights <- (snp_set$rda_model$CCA$eig/sum(snp_set$rda_model$CCA$eig))[1:K]list_AI <-lapply(list(clim_ref,clim_pred), function(df){lapply(1:K, function(i){# Below we calculate the adaptive index for the axis i# For that, we multiply the value of the variables at the location of the population# by the loadings of the variables along the axis iai <- df %>% dplyr::select(any_of(clim_var)) %>%apply(1, function(x) sum(x * snp_set$rda_model$CCA$biplot[,i]))if(!is.null(weights)){ai <- ai * weights[i]} }) %>%setNames(paste0("RDA", 1:length(.))) %>%as.data.frame()})if(CG==F){eucloffset <-unlist(lapply(1:nrow(list_AI[[1]]), function(x) dist(rbind(list_AI[[1]][x,], list_AI[[2]][x,]), method ="euclidean")))} elseif(CG==T){eucloffset <-lapply(1:nrow(list_AI[[2]]), function(x) as.matrix(dist(rbind(list_AI[[2]][x,],list_AI[[1]]), method ="euclidean"))[2:(nrow(list_AI[[1]])+1),1]) %>%setNames(clim_pred[,1] %>%pull()) %>%as.data.frame() %>%cbind(clim_ref[,1]) }return(eucloffset)}
# Function to make the genomic offset mapssource(here("scripts/functions/make_go_map.R"))# Population coordinatespop_coord <-readRDS(here(paste0("data/ClimaticData/MaritimePinePops/ClimatePopulationLocationPointEstimates_ReferencePeriods_noADJ.rds")))[[1]]$ref_means %>% dplyr::select(pop,longitude,latitude)# Find minimum and maximum values of genomic offset for the mapsgo_limits <-lapply(snp_sets, function(x) {lapply(names(list_dist_env), function(gcm){x$go[[gcm]]}) %>%unlist()}) %>%unlist() %>%range()# The minimum GO value is very very small, almost zero, so we fix it to zero.go_limits[[1]] <-0# Generate the maps for each set of SNPs and each GCMlapply(snp_sets, function(x) {go_maps <-lapply(names(list_dist_env), function(gcm){make_go_map(dfcoord=pop_coord,snp_set = x,gcm=gcm,ggtitle=gcm,go_limits = go_limits,point_size =3)})legend_maps <-get_legend(go_maps[[1]])go_maps <-lapply(go_maps, function(y) y +theme(legend.position ="none"))go_maps$legend_maps <- legend_mapsgo_maps <-plot_grid(plotlist=go_maps)# =====================================================# We save the figures for the Supplementary Information# =====================================================if(x$set_code=="cand"){ ggsave(here("figs/RDA/GOMaps_PopLocations_CandidateSNPs_SI.pdf"), go_maps, width=10,height=6, device="pdf")} elseif(x$set_code=="control"){ggsave(here("figs/RDA/GOMaps_PopLocations_ControlSNPs_SI.pdf"), go_maps, width=10,height=6, device="pdf")}# =========# Add title# =========title <-ggdraw() +draw_label( x$set_name,fontface ='bold',x =0,hjust =0 ) +theme(plot.margin =margin(0, 0, 0, 7))# merge title and plotsplot_grid( title, go_maps,ncol =1,rel_heights =c(0.1, 1)) })
For each GCM, we attribute the value 1 to the top five populations with the highest genomic offset and we attribute the value 0 to the other populations. We then count the number of 1 for each population, which gives the table and map below:
Code
source(here("scripts/functions/make_high_go_pop_maps.R"))high_go_pops <-make_high_go_pop_maps(pop_coord=pop_coord,list_go = snp_sets$cand$go,ggtitle="RDA",nb_id_pop =5) # number of selected populationssaveRDS(high_go_pops, file =here("outputs/RDA/high_go_pops.rds"))high_go_pops[[1]] %>% kable_mydf
pop
longitude
latitude
GFDL-ESM4
IPSL-CM6A-LR
MPI-ESM1-2-HR
MRI-ESM2-0
UKESM1-0-LL
sum_go
SEG
-8.45
42.82
0
0
0
0
0
0
CAD
-6.42
43.54
0
0
0
0
0
0
ARM
-6.46
43.30
0
0
0
0
0
0
SAC
-8.36
42.12
0
0
0
0
0
0
LAM
-6.22
43.56
0
0
0
0
0
0
SIE
-6.49
43.53
0
0
0
0
0
0
PUE
-6.63
43.55
0
0
0
0
0
0
ALT
-6.49
43.28
0
0
0
0
0
0
CAS
-6.98
43.50
0
0
0
0
0
0
OLB
-0.62
40.17
0
0
0
0
0
0
MIM
-1.30
44.13
0
0
0
0
0
0
PET
-1.30
44.06
0
0
0
0
0
0
PLE
-2.34
47.78
0
0
0
0
0
0
HOU
-1.15
45.18
0
0
0
0
0
0
STJ
-2.03
46.76
0
0
0
0
0
0
SAL
-3.06
41.83
0
0
0
0
0
0
VER
-1.09
45.55
0
0
0
0
0
0
TAM
-5.02
33.60
0
0
0
0
0
0
OLO
-1.83
46.57
0
0
0
0
0
0
BAY
-2.88
41.52
0
0
0
0
0
0
PIE
9.04
41.97
0
0
0
0
0
0
CAR
-4.28
41.17
0
0
0
0
0
0
LEI
-8.96
39.78
0
0
0
1
0
1
CEN
-4.49
40.28
0
0
1
0
0
1
VAL
-4.31
40.52
0
0
1
0
0
1
BON
-1.66
39.99
0
0
0
0
1
1
ORI
-2.35
37.53
0
0
0
0
1
1
QUA
-0.36
38.97
0
0
0
1
1
2
ARN
-5.12
40.19
0
0
1
1
0
2
COC
-4.50
41.25
1
1
0
0
0
2
CUE
-4.48
41.37
1
1
0
0
0
2
PIA
9.46
42.02
1
1
0
0
0
2
MAD
-5.23
35.18
1
1
1
1
1
5
COM
-3.95
36.83
1
1
1
1
1
5
Code
high_go_pops[[2]]
Comparing GO predictions
We look at the correlation across the different genomic offset predictions at the location of the populations, i.e. those based on all SNPs and those based on sets of candidates or control SNPs.
# Function to calculate the RDA genetic offset using rasters# ==========================================================# Arguments# =========# snp_set = list with at least the RDA models # clim_var = selected climatic variables# K = number of RDA axes used to calculate the genomic offset# gcm = Global Climate Model of the raster with future climatic data# clim_ref_adj = TRUE or FALSE, specify whether the point estimate climatic data used to scale the rasters should be adjusted or not for elevation# ref_period = the reference period used to calculate the adaptive index, can be 1901-1950 or 1960-1991# method = `loadings` or `predict` pred_GO_rasters <-function(snp_set, clim_var, K, gcm,range_buffer =shapefile(here('data/Mapping/PinpinDistriEUforgen_NFIplotsBuffer10km.shp')),method ="loadings",clim_ref_adj =FALSE,ref_period ="ref_1901_1950"){# Load point estimate climatic data of the reference periodif(clim_ref_adj==TRUE){adj <-"ADJ"} else {adj <-"noADJ"} clim_ref_pe <-readRDS(here(paste0("data/ClimaticData/MaritimePinePops/ClimatePopulationLocationPointEstimates_ReferencePeriods_",adj,".rds")))[[ref_period]]# Extract scaling parameters, i.e. mean and variance scale_params <-lapply(clim_var, function(x){ vec_var <- clim_ref_pe$ref_means[,x] %>%pull()list(mean =mean(vec_var), sd =sd(vec_var))}) %>%setNames(clim_var)# Rasters with climatic data for the reference periodpath <-here(paste0("data/ClimaticData/ClimateDTRasters/1km_",clim_ref_pe$range[[1]],"-",clim_ref_pe$range[[2]],"_Extent-JulietteA/"))tif_paths <-lapply(clim_var, function(x) paste0(path,"/",x,".tif"))ref_rasts <- raster::stack(tif_paths)# Raster with future climatic datapath <-here(paste0("data/ClimaticData/ClimateDTRasters/1km_",gcm,"_2041-2070_ssp370_Extent-JulietteA/"))tif_paths <-lapply(clim_var, function(x) paste0(path,"/",x,".tif"))fut_rasts <- raster::stack(tif_paths)# checking that the CRS is the same for the buffer and the rastersif(identical(crs(range_buffer),crs(ref_rasts))==FALSE){stop(paste0("CRS of the range buffer is not the same as the CRS of the reference period raster."))} if(identical(crs(range_buffer),crs(fut_rasts))==FALSE){stop(paste0("CRS of the range buffer is not the same as the CRS of future climate rasters."))} # Mask with the range if suppliedif(!is.null(range_buffer)){ ref_rasts <-mask(ref_rasts, range_buffer) fut_rasts <-mask(fut_rasts, range_buffer)}# extract coordinates and climatic values, and mean-center the climatic dataref_vals <-as.data.frame(rasterToPoints(ref_rasts))for(i in clim_var){ref_vals[,i] <- (ref_vals[,i] - scale_params[[i]]$mean) / scale_params[[i]]$sd}fut_vals <-as.data.frame(rasterToPoints(fut_rasts))for(i in clim_var){fut_vals[,i] <- (fut_vals[,i] - scale_params[[i]]$mean) / scale_params[[i]]$sd}# Predict genomic offset for each pixel based on the loadings of the climatic variablesif(method =="loadings"){ ref_proj <-list() fut_proj <-list() go_proj <-list()# Projection for each RDA axisfor(i in1:K){# Calculate adaptive index for the reference period ref_rast <- ref_rasts[[1]] ref_rast[!is.na(ref_rast)] <-as.vector(apply(ref_vals[,clim_var], 1, function(x) sum( x * snp_set$rda_model$CCA$biplot[,i])))names(ref_rast) <-paste0("RDA_ref_", as.character(i)) ref_proj[[i]] <- ref_rastnames(ref_proj)[i] <-paste0("RDA", as.character(i))# Calculate adaptive index under future climates fut_rast <- fut_rasts[[1]] fut_rast[!is.na(fut_rast)] <-as.vector(apply(fut_vals[,clim_var], 1, function(x) sum( x * snp_set$rda_model$CCA$biplot[,i])))names(fut_rast) <-paste0("RDA_fut_", as.character(i)) fut_proj[[i]] <- fut_rastnames(fut_proj)[i] <-paste0("RDA", as.character(i))# Calculate genetic offset based on a single RDA axis go_proj[[i]] <-abs(ref_proj[[i]] - fut_proj[[i]])names(go_proj)[i] <-paste0("RDA", as.character(i)) }}# Predict genomic offset for each pixel based on predict.RDAif(method =="predict"){# Prediction with the RDA model under reference-period climates and future climates ref_pred <-predict(snp_set$rda_model, ref_vals[,-c(1,2)], type ="lc") fut_pred <-predict(snp_set$rda_model, fut_vals[,-c(1,2)], type ="lc") ref_proj <-list() fut_proj <-list() go_proj <-list()for(i in1:K){# Calculate adaptive index for the reference period ref_rast <-rasterFromXYZ(data.frame(ref_vals[,c(1,2)], Z =as.vector(ref_pred[,i])), crs =crs(ref_rasts))names(ref_rast) <-paste0("RDA_ref_", as.character(i)) ref_proj[[i]] <- ref_rastnames(ref_proj)[i] <-paste0("RDA", as.character(i))# Calculate adaptive index under future climates fut_rast <-rasterFromXYZ(data.frame(ref_vals[,c(1,2)], Z =as.vector(fut_pred[,i])), crs =crs(ref_rasts))names(fut_rast) <-paste0("RDA_fut_", as.character(i)) fut_proj[[i]] <- fut_rastnames(fut_proj)[i] <-paste0("RDA", as.character(i))# Calculate genetic offset based on a single RDA axis go_proj[[i]] <-abs(ref_proj[[i]] - fut_proj[[i]])names(go_proj)[i] <-paste0("RDA", as.character(i)) } }# Weights based on axis eigen values weights <- snp_set$rda_model$CCA$eig/sum(snp_set$rda_model$CCA$eig)# Weight the current and future adaptive indices based on the eigen values of the associated axes ref_proj_weighted <-do.call(cbind, lapply(1:K, function(x) rasterToPoints(ref_proj[[x]])[,-c(1,2)])) ref_proj_weighted <-as.data.frame(do.call(cbind, lapply(1:K, function(x) ref_proj_weighted[,x]*weights[x]))) fut_proj_weighted <-do.call(cbind, lapply(1:K, function(x) rasterToPoints(fut_proj[[x]])[,-c(1,2)])) fut_proj_weighted <-as.data.frame(do.call(cbind, lapply(1:K, function(x) fut_proj_weighted[,x]*weights[x])))# Predict a global genetic offset, incorporating the K first axes weighted by their eigen values go_global_proj <- go_proj[[1]] go_global_proj[!is.na(go_global_proj)] <-unlist(lapply(1:nrow(ref_proj_weighted), function(x) dist(rbind(ref_proj_weighted[x,], fut_proj_weighted[x,]), method ="euclidean")))names(go_global_proj) <-"global_offset"# Return projections for current and future climates for each RDA axis, prediction of genetic offset for each RDA axis and a global genetic offset return(list(ref_proj = ref_proj, fut_proj = fut_proj, go_proj = go_proj, go_global_proj = go_global_proj, weights = weights[1:K]))}
Code
# For each snp set and each GCM, we generate a list of rasters with:# the projection of AI for reference_period and future climates# genomic offset predictions for each RDA axis # genomic offset predictions incorporating the K first axes weighted by their eigen valuesgo_rasters <- snp_sets %>%lapply(function(snp_set){lapply(names(clim_dfs$clim_pred), function(gcm){pred_GO_rasters(snp_set,clim_var,K=2,gcm)}) %>%setNames(names(clim_dfs$clim_pred)) })go_rasters %>%saveRDS(file=here("outputs/RDA/go_rasters.rds"))
Code
go_rasters <-readRDS(file=here("outputs/RDA/go_rasters.rds"))# Map options# ===========point_size=2x_limits =c(-10, 15)y_limits =c(31, 52)# get maximum and minimum genomic offset values across the different GO predictions (for each snp set and each GCM)# so that comparing the maps will be easiergo_limits <-names(snp_sets) %>%lapply(function(x){lapply(names(clim_dfs$clim_pred), function(gcm){rasterToPoints(go_rasters[[x]][[gcm]]$go_global_proj) %>%as_tibble() %>%pull(global_offset)}) %>%unlist()}) %>%unlist() %>%range()# The minimum GO value is very very small, almost zero, so we fix it to zero.go_limits[[1]] <-0# Mapping# =======go_maps <-names(snp_sets) %>%lapply(function(x){go_maps <-lapply(names(clim_dfs$clim_pred), function(gcm){# extract genomic offset values go_dfs <-rasterToPoints(go_rasters[[x]][[gcm]]$go_global_proj) %>%as_tibble()ggplot(data = go_dfs) +geom_sf(data = world, fill="gray98") +scale_x_continuous(limits = x_limits) +scale_y_continuous(limits = y_limits) +geom_raster(aes(x = x, y = y, fill = global_offset), alpha =1) +scale_fill_gradient2(low="blue", mid="yellow", high="red", midpoint=(max(go_limits[[2]])-min(go_limits[[1]]))/2,limits=go_limits,name ="Genomic offset") +xlab("Longitude") +ylab("Latitude") +ggtitle(gcm) +theme_bw(base_size =11) +theme(panel.grid =element_blank(), plot.background =element_blank(), panel.background =element_blank(), #legend.position = c(0.85,0.15),legend.box.background =element_rect(colour ="gray80"),strip.text =element_text(size=11))})legend_maps <-get_legend(go_maps[[1]])go_maps <-lapply(go_maps, function(y) y +theme(legend.position ="none"))go_maps$legend_maps <- legend_mapsgo_maps <-plot_grid(plotlist=go_maps)# We save the maps for the Supplementary Information# ==================================================if(x=="cand"){ ggsave(here("figs/RDA/GOMap_CandidateSNPs_SI.pdf"), go_maps, width=11,height=7, device="pdf")} elseif(x=="control"){ggsave(here("figs/RDA/GOMap_ControlSNPs_SI.pdf"), go_maps, width=11,height=7, device="pdf")}# Add title# =========title <-ggdraw() +draw_label( snp_sets[[x]]$set_name,fontface ='bold',x =0,hjust =0) +theme(plot.margin =margin(0, 0, 0, 7))# merge title and plotsplot_grid(title, go_maps, ncol =1,rel_heights =c(0.1, 1))})# We save the maps for the Github repository# ==========================================pdf(here("figs/RDA/GOmaps.pdf"), width=14,height=9.5)lapply(go_maps, function(x) x)dev.off()# Show GO maps# ============lapply(go_maps, function(x) x)
Corr btw predictions with rasters or spatial points
We check that the genomic offset predictions we obtained with the spatial points are highly correlated (ie similar) to those obtained with the rasters.
Comment When we extract the genomic offset values from the rasters at the location of the populations, we have missing values for some populations. Indeed, some populations are not included in the buffer of the species range that I used, based on the EUFORGEN distribution and 10km around the NFI plots. We may want to modify the buffer of the species range to include those populations!
# Load the climatic data of the NFI plots.nfi_clim <-readRDS(here("data/ClimaticData/NFIplots/NFIclimate.rds"))# Keep only the climatic variables of interest and scale the climatic datanfi_dfs <-generate_scaled_clim_datasets(clim_var, clim_ref = nfi_clim$clim_ref, clim_pred = nfi_clim$clim_survey)# calculate the genomic offset for the NFI plotssnp_sets <-lapply(snp_sets, function(x){x$go_nfi <-pred_GO_spatial_points(snp_set = x,K =2, # why K=2??clim_var=clim_var,clim_ref = nfi_dfs$clim_ref,clim_pred = nfi_dfs$clim_pred,weights = T)return(x)})# checking missing data# lapply(snp_sets, function(x) sum(is.na(x$go_nfi)))# Find minimum and maximum values of genomic offset for the mapsgo_limits <-lapply(snp_sets, function(snp_set) snp_set$go_nfi) %>%unlist() %>%range()# The minimum GO value is very very small, almost zero, so we fix it to zero.go_limits[[1]] <-0# map genomic offset predictions in the NFI plots lapply(snp_sets, function(x) {p <-make_go_map(dfcoord=readRDS(here("data/ClimaticData/NFIplots/NFIclimate.rds"))[[1]] %>% dplyr::select(contains("ude")), snp_set = x,type="NFI",point_size =0.5,go_limits = go_limits,legend_position =c(0.85,0.15),legend_box_background ="gray80",y_limits =c(35, 51))# Figure for the SI# =================if(x$set_code=="cand"){ p_SI <- p +theme(plot.title =element_blank())ggsave(here("figs/RDA/NFI_GOmap_CandidateSNPs_SI.pdf"), p_SI, width=8,height=8, device="pdf") } elseif(x$set_code=="control"){ p_SI <- p +theme(plot.title =element_blank())ggsave(here("figs/RDA/NFI_GOmap_ControlSNPs_SI.pdf"), p_SI, width=8,height=8, device="pdf") }# Show maps in the Quarto document# ================================p})
We look at the correlation across the different genomic offset predictions in the NFI plots, i.e. those based on all SNPs and those based on sets of candidates or control SNPs.
cg_clim <-readRDS(here("data/ClimaticData/CommonGardens/ClimateCG.rds")) %>% dplyr::select(cg,any_of(clim_var))cg_coord <-readRDS(here("data/ClimaticData/CommonGardens/ClimateCG.rds")) %>% dplyr::select(cg,contains("ude"))cg_names <-unique(cg_coord$cg)# Generate scaled climatic datasets with climatic data at the location of the populations and at the location of the common gardenscg_dfs <-generate_scaled_clim_datasets(clim_var, clim_pred = cg_clim)# Predict genomic offset of each population when transplanted in the climate of the common gardenssnp_sets <-lapply(snp_sets, function(x){x$go_cg <-pred_GO_spatial_points(snp_set = x,K =2, # why K=2??clim_var = clim_var,clim_ref = cg_dfs$clim_ref,clim_pred = cg_dfs$clim_pred,weights = T,CG=T) %>%as_tibble()return(x)})# Map genomic offset predictions at the locations of the populationsgo_maps_cg <-lapply(cg_names, function(cg_name){# Find minimum and maximum values of genomic offset for the mapsgo_limits <-lapply(snp_sets, function(snp_set) snp_set$go_cg[[cg_name]]) %>%unlist() %>%range()go_limits[[1]] <-0p <-lapply(snp_sets, function(x) make_go_map(dfcoord=pop_coord, snp_set = x,point_size =3,type="CG",go_limits = go_limits,cg_name=cg_name,cg_coord=cg_coord))plot_grid(p[[1]],p[[2]],p[[3]],nrow=1)}) %>%setNames(cg_names)pdf(here("figs/RDA/GOmaps_CGs.pdf"), width=15,height=6)lapply(go_maps_cg, function(x) x)dev.off()# show mapslapply(go_maps_cg, function(x) x)
We look at the correlation across the different genomic offset predictions in the common gardens, i.e. those based on all SNPs and those based on sets of candidates or control SNPs.
Capblancq, Thibaut, and Brenna R Forester. 2021. “Redundancy Analysis (RDA): A Swiss Army Knife for Landscape Genomics.”Methods in Ecology and Evolution.
Capblancq, Thibaut, Susanne Lachmuth, Matthew C Fitzpatrick, and Stephen R Keller. 2023. “From Common Gardens to Candidate Genes: Exploring Local Adaptation to Climate in Red Spruce.”New Phytologist 237 (5): 1590–1605.